#!/usr/bin/perl
# 2014/01/30 大幅な改訂。
# 2011/02/05 automation program for ChIP-seq
# ver2 controlとsampleで数が違う場合に、sampleに合わせる。足りない場合には追加、多い場合には減らす
# ver3 controlとsampleで数が違う場合に、少ない方に合わせる。uniqueなタグが同じになるようにsamplingする。
# ver5 re-samplingのステップをとり除く
# ver5.1 Use MAQ instead of Bowtie (but convert to Bowtie format)

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use FindBin;
use Carp qw(croak);
$| = 0;


#---------------------------------------
# このプログラムの実行前にしておくこと
#---------------------------------------
# 1. $ perl Batch_FastQ2Bam.pl -i [input fastq file]
#     1.1 BWAでアラインメント
#     1.2 samファイルをbamファイルに変換
#     1.3 bamファイルをソート
#     1.4 indexを作成
# 2. $ export PATH=/cscb/cscbsoftware/centos5-x86_64/local/homer/bin/:$PATH
#    $ batchMakeTagDirectories.pl key.txt -cpu 12 -keepAll 2>&1 | tee MakeTagDirectory.log
#    (key.txt の中身は、"tagDirectoryName <tab> path-to-alignment-file"という形式の行を必要な分だけ列挙したもの)
#     2.1 tsvファイルを作成する
#     2.2 ピークシフトを予測する

#---------------------------------------
# parameter
#---------------------------------------
my %param = (
	6	=> 1,
	8	=> 1,
);
my $program_DIR = $FindBin::Bin;


if( !(exists $param{scalar(@ARGV)}) or $ARGV[0] eq '--help'){
	my $help =<<"EOF";
Usage : $0 -s [treatment directory] -c [control directory] -o [output directory] -j [jump position]
RATIO: background calculation
SIGNAL_MAP: create signal map
SCORE: score calculation
MERGE: merge peaks
EOF
	die "$help";
}

### get option
my %opt;
getopts("s:c:j:o:", \%opt);
my $DIR_treatment = $opt{s};
$DIR_treatment =~ s|/$||;
my $DIR_control = $opt{c};
$DIR_control =~ s|/$||;
my $DIR_output = $opt{o};
$DIR_output =~ s|/$||;
my $label = $opt{j};

#---------------------------------------
# parameter
#---------------------------------------
my $FILE_window_reads = $DIR_output . '/' . '1kb_window_comparison.txt';
my $FILE_back_ratio = $DIR_output . '/' . 'background_ratio.txt';
my $FILE_score = $DIR_output . '/' . 'score.txt';
my $FILE_score_compressed = $DIR_output . '/' . 'score.txt.gz';
my $FILE_significant_region = $DIR_output . '/' . 'significant_region.txt';
my $FILE_mergethreshold = $DIR_output . '/' . 'merge_threshold.txt';
my $FILE_window_plot = $DIR_output . '/' . 'window_plot.jpg';
my $FILE_average_treatment_score_of_chrI_chrII = $DIR_output . '/' . 'average_treatment_score_of_chrI_chrII.txt';

printf "[%s] start automation process\n", &time();

### jump to some position
if(exists $opt{j}){
	goto $label;
}

#---------------------------------------
# background ratio calculation
#---------------------------------------
&cmd("background_ratio.pl -s $DIR_treatment -c $DIR_control -o $FILE_back_ratio", "ratio calculation");


#---------------------------------------
# make signal map and wig file (100bp widow, 30bp sliding)
#---------------------------------------
SIGNAL_MAP:
{
	my $f1 = $DIR_treatment . '/' . 'I.sgr';
	unless(-e $f1){
		&cmd("signalmap.pl -d $DIR_treatment", "signal map for treatment");
	}
	my $f2 = $DIR_control . '/' . 'I.sgr';
	unless(-e $f2){
		&cmd("signalmap.pl -d $DIR_control", "signal map for control");
	}
}

#---------------------------------------
# calculate score & calculate merge threshold
#---------------------------------------
SCORE:
my $background_ratio = `cat $FILE_back_ratio`;
$background_ratio =~ s/\r?\n//;
&Rcmd("score_calculation.R $DIR_treatment $DIR_control $background_ratio $FILE_score $FILE_mergethreshold $FILE_average_treatment_score_of_chrI_chrII", "calc score & calc merge threshold");
system("gzip $FILE_score") == 0 or die "compress failed: $!";

#---------------------------------------
# merge peaks
#---------------------------------------
MERGE:
my $threshold = `cat $FILE_mergethreshold`;
$threshold =~ s/\r?\n//;
&cmd("mergePeaks.pl -i $FILE_score_compressed -t $threshold -o $FILE_significant_region", "merge peaks");

#---------------------------------------
# calculate P-value & Q-value
#---------------------------------------
&Rcmd("Statistics.R $FILE_significant_region", "P-value & Q-value calculation");



sub cmd{
	my ($command, $message) = @_;
	my $time = &time();
	print "[$time] $message ...";
	system("perl $program_DIR/$command") == 0 or die "$command faild: $!";
	print "OK\n";
}

sub Rcmd{
	my ($command, $message) = @_;
	my $time = &time();
	print "[$time] $message ...";
	system("Rscript --vanilla --slave $program_DIR/$command") == 0 or die "$command failed: $!";
	print "OK\n";
}


sub time{
	my ($sec, $min, $hour, $mday, $mon, $year) = localtime(time);
	my $time = sprintf ("%d/%0d %02d:%02d",$mon+1, $mday, $hour, $min);
	return $time;
}
